import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings

warnings.filterwarnings('ignore')

# ====================================================
# 数据准备
# ====================================================

# 读取历史销售数据
sales_data = pd.read_excel(r"E:\M\MathModel\23国赛C题\ds\23年6月24到30.xlsx")
# 读取回归分析结果
regression_results = pd.read_excel(r"E:\M\MathModel\23国赛C题\ds\销售分析图表\回归分析结果.xlsx")
# 读取补货量预测
replenishment_pred = pd.read_excel(r"E:\M\MathModel\23国赛C题\ds\7月1日补货量预测\7月1日单品补货量预测.xlsx")


# 数据预处理
def preprocess_data(sales_data, regression_results, replenishment_pred):
    """整合并预处理数据"""
    # 获取批发价（进价）
    wholesale_prices = sales_data.groupby('单品名称')['批发价'].mean()

    # 合并数据
    merged_data = pd.merge(replenishment_pred, wholesale_prices, on='单品名称', how='left')
    merged_data = pd.merge(merged_data, regression_results, on='单品名称', how='left')

    # 重命名列
    merged_data.rename(columns={
        '预测销量(kg)': 'predicted_sales',
        '损耗率(%)': 'loss_rate',
        '补货量(kg)': 'replenishment',
        '批发价': 'wholesale_price',
        '斜率': 'slope',
        '截距': 'intercept'
    }, inplace=True)

    # 转换损耗率为小数形式
    merged_data['loss_rate'] = merged_data['loss_rate'] / 100.0

    # 处理斜率为0的情况（价格不变）
    merged_data['slope'] = merged_data['slope'].replace(0, 1e-10)

    # 确保所有必要列都有值
    merged_data.fillna({
        'wholesale_price': merged_data['wholesale_price'].mean(),
        'slope': merged_data['slope'].mean(),
        'intercept': merged_data['intercept'].mean()
    }, inplace=True)

    # 添加验证：检查价格是否合理
    merged_data['min_price'] = merged_data.apply(
        lambda row: min(row['intercept'], row['intercept'] + row['slope'] * row['replenishment']),
        axis=1
    )
    merged_data['max_price'] = merged_data.apply(
        lambda row: max(row['intercept'], row['intercept'] + row['slope'] * row['replenishment']),
        axis=1
    )

    # 打印有问题的商品
    problem_products = merged_data[
        (merged_data['min_price'] < 0) |
        (merged_data['max_price'] > 100) |
        (merged_data['slope'].abs() > 10)
        ]

    if not problem_products.empty:
        print("\n警告：以下商品的价格-销售量关系可能存在问题:")
        print(problem_products[['单品名称', 'slope', 'intercept', 'min_price', 'max_price']])

    return merged_data


# 获取处理后的数据
product_data = preprocess_data(sales_data, regression_results, replenishment_pred)
print("数据准备完成，共处理{}种商品".format(len(product_data)))
print(product_data[['单品名称', 'slope', 'intercept', 'wholesale_price', 'replenishment']].head())


# ====================================================
# 非线性规划模型
# ====================================================

class ProfitOptimizationModel:
    def __init__(self, product_data):
        """初始化优化模型"""
        self.products = product_data
        self.num_products = len(product_data)

    def calculate_profit(self, X):
        """
        计算总利润
        X: 各商品的销售量数组 (1D numpy array)
        """
        total_profit = 0.0
        for i, row in self.products.iterrows():
            # 获取商品参数
            a = row['intercept']
            b = row['slope']
            T = row['replenishment']
            J = row['wholesale_price']
            K = row['loss_rate']

            # 计算售价
            P = a + b * X[i]

            # 计算成本变动总额
            B = K * P * T

            # 计算单品利润并累加
            profit_i = P * X[i] - (B + T * J)
            total_profit += profit_i

        return total_profit

    def check_constraints(self, X):
        """
        检查约束条件
        返回: (是否满足所有约束, 违反约束的商品索引)
        """
        violations = []
        for i, row in self.products.iterrows():
            # 获取商品参数
            a = row['intercept']
            b = row['slope']
            T = row['replenishment']
            J = row['wholesale_price']

            # 计算售价
            P = a + b * X[i]

            # 检查约束条件
            if X[i] > T:  # 销售量不超过补货量
                violations.append(i)
            if P <= J:  # 售价应大于进价
                violations.append(i)
            if T < 2.5:  # 最小陈列量
                violations.append(i)

        return (len(violations) == 0, violations)


# 初始化模型
optimization_model = ProfitOptimizationModel(product_data)


# ====================================================
# 改进的粒子群算法 (PSO) - 修复边界计算问题
# ====================================================

class ImprovedPSO:
    def __init__(self, model, num_particles=30, max_iter=100,
                 w_max=0.9, w_min=0.4, c1=2.0, c2=2.0, sigma=0.1):
        """初始化改进的PSO算法"""
        self.model = model
        self.num_particles = num_particles
        self.max_iter = max_iter
        self.w_max = w_max
        self.w_min = w_min
        self.c1 = c1
        self.c2 = c2
        self.sigma = sigma
        self.num_dim = model.num_products

        # 合理的边界计算
        self.lower_bounds = np.zeros(self.num_dim)  # 最小销售量设为0
        self.upper_bounds = np.zeros(self.num_dim)  # 初始化上界

        # 计算每个商品的上界
        for i, row in enumerate(model.products.iterrows()):
            _, product = row
            slope = product['slope']
            intercept = product['intercept']
            wholesale_price = product['wholesale_price']
            replenishment = product['replenishment']

            # 对于斜率极小的情况，限制在合理范围内
            if abs(slope) < 1e-5:
                self.lower_bounds[i] = 0
                self.upper_bounds[i] = replenishment
                continue

            # 计算理论边界
            if slope > 0:
                # 价格随销售量增加而上升
                critical_x = (wholesale_price - intercept) / slope
                self.lower_bounds[i] = max(0, min(critical_x, replenishment))
                self.upper_bounds[i] = replenishment
            else:
                # 价格随销售量增加而下降
                critical_x = (intercept - wholesale_price) / (-slope)
                self.upper_bounds[i] = min(replenishment, critical_x)
                self.lower_bounds[i] = 0

            # 确保边界合理
            self.lower_bounds[i] = max(0, self.lower_bounds[i])
            self.upper_bounds[i] = min(replenishment, self.upper_bounds[i])

            # 防止边界无效
            if self.upper_bounds[i] < self.lower_bounds[i]:
                self.upper_bounds[i] = self.lower_bounds[i] = replenishment * 0.5

        # 打印边界信息
        print("\n变量边界信息:")
        for i, (low, high) in enumerate(zip(self.lower_bounds, self.upper_bounds)):
            product = model.products.iloc[i]
            print(f"商品 {i + 1} ({product['单品名称']}):")
            print(f"  下界={low:.2f}, 上界={high:.2f} (补货量={product['replenishment']:.2f})")
            print(
                f"  斜率={product['slope']:.6f}, 截距={product['intercept']:.2f}, 进价={product['wholesale_price']:.2f}")

        # 初始化粒子群
        self.positions = np.zeros((num_particles, self.num_dim))
        self.velocities = np.zeros((num_particles, self.num_dim))
        self.personal_best_positions = np.zeros((num_particles, self.num_dim))
        self.personal_best_scores = np.full(num_particles, -np.inf)

        # 全局最优
        self.global_best_position = np.zeros(self.num_dim)
        self.global_best_score = -np.inf

        # 记录收敛过程
        self.convergence = []

        # 初始化粒子
        self.initialize_particles()

    def initialize_particles(self):
        """初始化粒子位置和速度"""
        for i in range(self.num_particles):
            # 在边界内随机初始化位置
            self.positions[i] = np.array([
                np.random.uniform(low, high)
                for low, high in zip(self.lower_bounds, self.upper_bounds)
            ])

            # 初始化速度
            self.velocities[i] = np.random.uniform(
                -0.1 * (self.upper_bounds - self.lower_bounds),
                0.1 * (self.upper_bounds - self.lower_bounds)
            )

            # 计算初始适应度
            score = self.model.calculate_profit(self.positions[i])

            # 更新个体最优
            self.personal_best_positions[i] = self.positions[i].copy()
            self.personal_best_scores[i] = score

            # 更新全局最优
            if score > self.global_best_score:
                self.global_best_position = self.positions[i].copy()
                self.global_best_score = score

        print(f"初始全局最优利润: {self.global_best_score:.2f}")

    def update_inertia_weight(self, iter):
        """动态更新惯性权重"""
        return self.w_max - (self.w_max - self.w_min) * (iter / self.max_iter)

    def run(self):
        """执行PSO优化"""
        for iter in range(self.max_iter):
            # 更新惯性权重
            w = self.update_inertia_weight(iter)

            for i in range(self.num_particles):
                # 生成随机因子
                r1 = np.random.rand(self.num_dim)
                r2 = np.random.rand(self.num_dim)

                # 计算随机惯性权重
                rand_w = w + (w - self.w_min) * np.random.rand() + self.sigma * np.random.randn()

                # 更新速度
                cognitive = self.c1 * r1 * (self.personal_best_positions[i] - self.positions[i])
                social = self.c2 * r2 * (self.global_best_position - self.positions[i])
                self.velocities[i] = rand_w * self.velocities[i] + cognitive + social

                # 限制速度范围
                velocity_max = 0.2 * (self.upper_bounds - self.lower_bounds)
                self.velocities[i] = np.clip(self.velocities[i], -velocity_max, velocity_max)

                # 更新位置
                self.positions[i] += self.velocities[i]

                # 边界处理
                self.positions[i] = np.clip(self.positions[i], self.lower_bounds, self.upper_bounds)

                # 计算新位置的适应度
                current_score = self.model.calculate_profit(self.positions[i])

                # 更新个体最优
                if current_score > self.personal_best_scores[i]:
                    self.personal_best_positions[i] = self.positions[i].copy()
                    self.personal_best_scores[i] = current_score

                    # 更新全局最优
                    if current_score > self.global_best_score:
                        self.global_best_position = self.positions[i].copy()
                        self.global_best_score = current_score

            # 记录当前迭代的最优值
            self.convergence.append(self.global_best_score)

            # 打印进度
            if (iter + 1) % 10 == 0:
                print(f"迭代 {iter + 1}/{self.max_iter}, 当前最优利润: {self.global_best_score:.2f}")

        return self.global_best_position, self.global_best_score


# ====================================================
# 模型求解与结果分析
# ====================================================

# 创建并运行PSO优化器
print("\n开始粒子群优化...")
pso = ImprovedPSO(optimization_model, num_particles=30, max_iter=100)
best_X, best_profit = pso.run()
print(f"\n优化完成，最大利润: {best_profit:.2f}")

# 计算最优售价
best_prices = []
for i, row in product_data.iterrows():
    P = row['intercept'] + row['slope'] * best_X[i]
    best_prices.append(P)

# 创建结果DataFrame
results_df = pd.DataFrame({
    '单品名称': product_data['单品名称'],
    '补货量(kg)': product_data['replenishment'],
    '销售价格': best_prices,
    '预测销售量': best_X,
    '批发价': product_data['wholesale_price'],
    '损耗率': product_data['loss_rate'],
    '单品利润': [P * X - (K * P * T + T * J)
                 for P, X, T, J, K in zip(best_prices, best_X,
                                          product_data['replenishment'],
                                          product_data['wholesale_price'],
                                          product_data['loss_rate'])]
})

# 打印部分结果（前3个商品）
print("\n优化结果（前3个商品）:")
print(results_df.head(3))

# 保存完整结果
output_path = r"E:\M\MathModel\23国赛C题\ds\7月1日单品定价优化结果.xlsx"
results_df.to_excel(output_path, index=False)
print(f"\n完整结果已保存至: {output_path}")

# ====================================================
# 结果可视化
# ====================================================

# 收敛曲线
plt.figure(figsize=(10, 6))
plt.plot(pso.convergence, 'b-', linewidth=2)
plt.xlabel('迭代次数')
plt.ylabel('总利润')
plt.title('粒子群算法收敛过程')
plt.grid(True)
plt.savefig(r"E:\M\MathModel\23国赛C题\ds\pso_convergence.png")
plt.show()

# 价格与销售量关系示例（第一个商品）
if len(product_data) > 0:
    product0 = product_data.iloc[0]
    a0, b0 = product0['intercept'], product0['slope']
    T0 = product0['replenishment']
    J0 = product0['wholesale_price']
    X_range = np.linspace(0.1, T0, 100)
    P_range = a0 + b0 * X_range

    # 计算利润曲线
    profit_range = [P * X - (product0['loss_rate'] * P * T0 + T0 * J0)
                    for P, X in zip(P_range, X_range)]

    plt.figure(figsize=(12, 8))

    # 价格-销售量关系
    plt.subplot(2, 1, 1)
    plt.plot(X_range, P_range, 'r-', label='价格-销售量关系')
    plt.axvline(x=best_X[0], color='b', linestyle='--', label='最优销售量')
    plt.axhline(y=best_prices[0], color='g', linestyle='--', label='最优价格')
    plt.axhline(y=J0, color='y', linestyle='-.', label='批发价')
    plt.xlabel('销售量 (kg)')
    plt.ylabel('销售价格 (元/kg)')
    plt.title(f"{product0['单品名称']} 价格-销售量关系与最优解")
    plt.legend()
    plt.grid(True)

    # 利润-销售量关系
    plt.subplot(2, 1, 2)
    plt.plot(X_range, profit_range, 'b-', label='利润曲线')
    plt.axvline(x=best_X[0], color='r', linestyle='--', label='最优销售量')
    plt.xlabel('销售量 (kg)')
    plt.ylabel('利润 (元)')
    plt.title(f"{product0['单品名称']} 利润-销售量关系")
    plt.legend()
    plt.grid(True)

    plt.tight_layout()
    plt.savefig(r"E:\M\MathModel\23国赛C题\ds\price_sales_profit_relationship.png")
    plt.show()

# 各商品利润贡献
if len(results_df) > 0:
    plt.figure(figsize=(14, 10))
    # 按利润排序
    sorted_results = results_df.sort_values('单品利润', ascending=True)
    # 只显示前20个商品
    display_count = min(20, len(sorted_results))
    plt.barh(sorted_results['单品名称'][-display_count:],
             sorted_results['单品利润'][-display_count:],
             color='skyblue')
    plt.xlabel('利润贡献 (元)')
    plt.title('各商品利润贡献')
    plt.tight_layout()
    plt.savefig(r"E:\M\MathModel\23国赛C题\ds\profit_contribution.png")
    plt.show()

# 价格与批发价比较
if len(results_df) > 0:
    plt.figure(figsize=(14, 8))
    n = len(results_df)
    index = np.arange(n)
    bar_width = 0.35

    plt.bar(index, results_df['批发价'], bar_width, label='批发价')
    plt.bar(index + bar_width, results_df['销售价格'], bar_width, label='优化售价')

    plt.xlabel('商品')
    plt.ylabel('价格 (元/kg)')
    plt.title('优化售价与批发价比较')
    plt.xticks(index + bar_width / 2, results_df['单品名称'], rotation=90)
    plt.legend()
    plt.tight_layout()
    plt.savefig(r"E:\M\MathModel\23国赛C题\ds\price_comparison.png")
    plt.show()

print("\n所有分析完成！")